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Abstract. We give a fast, exact algorithm for solving Dirichlet problems with 
polynomial boundary functions on quadratic surfaces in R," such as ellipsoids, 
elliptic cylinders, and paraboloids. To produce this algorithm, first we show 
that every polynomial in R" can be uniquely written as the sum of a harmonic 
function and a polynomial multiple of a quadratic function, thus extending 
a theorem of Ernst Fischer. We then use this decomposition to reduce the 
Dirichlet problem to a manageable system of linear equations. The algorithm 
requires differentiation of the boundary function, but no integration. We also 
show that the polynomial solution produced by our algorithm is the unique 
polynomial solution, even on unbounded domains such as elliptic cylinders and 
paraboloids. 



1. Introduction 

In this paper we present a fast, exact algorithm for solving Dirichlet problems 
with polynomial boundary functions on a quadratic surface in R" (n > 2). To 
illustrate the kind of Dirichlet problem we study, fix 6 = (6i, . . . , 6,i) £ R". For 
X = (xi, . . . , Xn) G R", we will write 

\\hxf^hlxl + --- + hlxl. 

Suppose we are given a polynomial p on R" . We wish to find a harmonic polynomial 
that equals p on the quadratic surface {x G R" : \\bx\\^ — 1}. 

Even if all the bj are nonzero (so that our quadratic surface is bounded and is, 
in particular, an ellipsoid), computing a solution to this Dirichlet problem presents 
several difficulties. A standard means of expressing the solution to the Dirichlet 
problem for bounded domains involves the Green's function and integration. How- 
ever, the Green's function of an ellipsoid does not have a known formula allowing 
for exact computations. An alternative approach avoids integration by employ- 
ing a finite difference, finite element, or Galerkin-type scheme to approximate the 
solution, but again this procedure will not produce an exact solution. 

If 6i = 62 = ■ • • = = 1, then our quadratic surface is the unit sphere. In this 
case, a fast algorithm for finding exact solutions is presented in |^ . That algorithm 
involves differentiation but no integration. The basis of that algorithm is that any 
polynomial p of degree m on R" can be decomposed in the form 

p = h+{\\x\\^-l)f, 

where ft, is a harmonic polynomial of degree at most m and / is a polynomial of 
degree at most m — 2. Because ft is harmonic and equals p on the unit sphere. 
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it is the solution to our Dirichlet problem. The algorithm presented in ||^ shows 
how the polynomials h and / in the decomposition above can be computed via 
differentiation from p. Unfortunately these techniques work only on spheres and so 
do not provide an algorithm for nonspherical ellipsoids or other quadratic surfaces. 

In this paper we solve the Dirichlet problem discussed above, getting solutions for 
ellipsoids as well as for elliptic cylinders (for example, {x G R'^ : x\ + 2x\ = 1}) and 
paraboloids (for example, {x G R'^ : x^ = xi -\- x\}). We will begin by extending 
the decomposition above to a collection of quadratic surfaces. We then use this 
decomposition to produce a system of linear equations whose solution will give an 
exact solution to our Dirichlet problem. We will show how this system of linear 
equations has a structure allowing it to be reduced to smaller systems of linear 
equations, thus producing a fast algorithm. The algorithm requires differentiation 
of the boundary function, but no integration. 

Before we turn to these matters, we need to present some background, much of 
which appears in A multi-index is an n-tuple a = (ai, . . . , a„) of nonnegative 
integers. The order of a, denoted |a|, is defined by 

|a| = ai + • • • + a„. 

We let x°' denote the monomial xi"^ . . . a;„"" and denote the differential oper- 
ator _Di"^ . . . Dn"" , where Dj denotes differentiation with respect to Xj. If g is a 
polynomial on R" given by q{x) = CqO:", then q{D) is the differential operator 
defined by q{D) = CaD°'. A polynomial is called homogeneous of degree m if 
it is a linear combination of monomials of degree m. 

Ernst Fischer proved that given a homogeneous polynomial q on R", every ho- 
mogeneous polynomial p of degree m can be decomposed uniquely as p = h + qf, 
where ft is a homogeneous polynomial of degree m satisfying q{D)h = and / is a 
homogeneous polynomial of suitable degree. 

In Q the subject of more general decompositions is discussed. Given two poly- 
nomials g and q on R", the relevant question is whether an arbitrary polynomial 
p can be decomposed as p — h -\- qf, where h and / are polynomials, with h satis- 
fying g{D)h — 0. Shapiro refers to a pair (gr, q) with this property as a generalized 
Fischer pair. He asks: Which {g, q) form generalized Fischer pairs? Note that if 
g{x) = then g{D) is the Laplacian and so this decomposition would require 

h to be harmonic. We will provide examples of a robust class of quadratic polyno- 
mials q that form a generalized Fischer pair with g{x) = Hxp, and we give explicit 
examples of the decomposition via our algorithm. 

Some of the surfaces that we consider are unbounded (for example, the elliptic 
cylinders and paraboloids mentioned above). Thus unique solutions to Dirichlet 
problems on these surfaces, even in the class of polynomials, are neither automatic 
nor expected. For example, the set of harmonic polynomials that vanish on the 
hyperplane {x G R" : Xn = 0} is not trivial. However, we will show that for the 
quadratic surfaces we consider, polynomial solutions to the Dirichlet problem with 
polynomial boundary functions are unique. 

The paper is organized as follows: In Section ^ we begin by presenting Fischer's 
lemma and the corresponding decomposition theorem. We then extend these re- 
sults to cover a wider class of generalized Fischer pairs. These generalized Fischer 
pairs are then used to solve the Dirichlet problem. We prove that the polynomial 
produced by this technique is the unique polynomial solving the Dirichlet prob- 
lem, even when our quadratic surface is unbounded. In Section H we describe a 
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fast algorithm for computing the solution to the Dirichlet problem promised by the 
results in Section |^. In Section ^ we present some examples, computed using our 
algorithm, of solutions to Dirichlet problems on ellipsoids. The Appendix contains 
a differentiation formula needed by our algorithm. 

2. Fischer's Lemma and the Dirichlet Problem 

In this section, we state and prove our generalization of Fischer's results. Then 
we prove the decomposition theorem that solves the Dirichlet problem. Then we 
show that even when our quadratic surfaces are unbounded, the solution given by 
our decomposition theorem is the only polynomial solution. 

We begin by stating Fischer's results, which are nicely restated and proved in 
also see Q. Fix an integer n > 2. We will always use m to denote a nonnegative 
integer. Let Vm denote the vector space of polynomials (with real coefficients) of 
degree at most m on R". For convenience, we define Vk to be {0} for fc < 0. As 
usual, A denotes the Laplacian. 

Lemma 2.1 (Fischer's Lemma). Suppose b = (6i,...,6„), where each bj ^ 0. 
Define L : Vm Vm by 

L(/) = A((||6xf -1)/). 
Then L is a linear, degree-preserving, bijection of Vm onto itself. 

Fischer's Lemma leads to Fischer's Decomposition Theorem, which gives a solu- 
tion to the Dirichlet problem for ellipsoids. 

Theorem 2.2 (Fischer's Decomposition Theorem). Suppose b — (6i, . . . , 6„), where 
each bj ^ 0. Let p e Vm- Then there exists a unique harmonic polynomial h G Vm 
such that 

p^h+i\\bx\\^-l)f 

for some f eVm~2- 

Let E — {x € R" : < 1}. Then h is the unique continuous function on E 

that is harmonic on E and equals p on the ellipsoid dE. 

We will need the following generalization of Corollary 5.3 of ||^. 

Lemma 2.3. If b — (5i, . . . , 6„), where each bj ^ 0, then no nonzero polynomial 
multiple of WbxW^ is harmonic. 

Proof. Suppose, to the contrary, that p is a nonzero polynomial of degree m such 
that llbxipp is harmonic. Let E be as in the theorem above, which states that there 
exists a harmonic polynomial h G Vm such that h equals p on the ellipsoid dE. But 
||6a;|pp is also a harmonic polynomial that equals p on dE. Because ||6a;|pp has 
degree m + 2, the polynomials and h cannot be equal. But this contradicts 

the uniqueness of solutions to the Dirichlet problem on bounded domains. □ 

In order to generate new generalized Fischer pairs, let g be a quadratic (degree 2) 
polynomial on R". We want to look at the map on Vm defined by 

L{f) = ^{qf). 

Our goal is to identify choices of q for which L is a bijection of Vm onto Vm- This 
leads us to the definition of a nonhyperbolic quadratic. 
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Definition 2.4. A nonhyperbolic quadratic is a polynomial q on R" of the form 



where at least one bj ^ 0. 



Note that the next theorem, which gives the desired bijectivity, imphes that no 
nonzero polynomial multiple of a nonhyperbolic quadratic is harmonic, generalizing 



Lemma 2.3. This result does not hold for arbitrary quadratic polynomials. In fact, 
even for a nonharmonic quadratic polynomial, a nonzero polynomial multiple might 
be harmonic. For example, — harmonic, but {xl — 3x2)xi is harmonic. 

Tiieorem 2.5. Let q be a nonhyperbolic quadratic. Define L: Vm ^ "Pm by 

L{f ) - A(<z/). 

Then L is a linear bijection of Vm onto Vm ■ 

Proof. Clearly L is a linear map of Vm into Vm- Since Vmi is finite dimensional, we 
need only show that L is injcctive. So suppose that this is not true. Then there 
exists / e Vm, / / 0, such that L{f) = 0. Without loss of generality, we may 
suppose that / is of degree m (otherwise, replace m by a lower integer). We can 
write / = /,„ + fm-i, where fm is homogeneous of degree m and fm-i G Vm-i- 
Write q = q2 + qi, where q2 is homogeneous degree 2 (so q2{x) = X]j=i ^'j^^ 
notation above) and qi £ Vi. Because ^{qf) — 0, we know that 

^{q-ifm + q-ifm-i + qif) = 0. 

Thus 

A(g2/m) = 0, 

because the other terms have lower degrees. We will now show that this implies 
that fm = 0, which is a contradiction. 

Note that we have reduced our theorem to the case where q = q2. Reordering 
the variables (if some of the bj =0), we see that it suffices to prove our theorem in 
the case when q{x) = X]j^=i ^i'^ii where 1 < r < n and 6i, . . . , 6^ are all nonzero. 
To simplify notation, we will also replace /,„ in the previous paragraph with /. So 
again we have the assumption that A(qf) = and we want to prove that / = 0, 
but now we have a special form for q. 



If r = n, our desired conclusion that f — follows from Lemma 2.2. So suppose 
r < n. Let k denote the degree of / thought of as a polynomial in Xr+i, . . . ,a;„ 
(temporarily think of xi , . . . , a;,, as constants to define k) . Write 

(2.6) f = P + g, 

where p is the part of / that is homogeneous of degree k in the variables Xr+i, . . . , a;„ 
and g is the remaining part of /, consisting of lower degree terms in the variables 

Using the product formula for the Laplacian, which states that 

(2.7) A(gfp) =pAgf + gfAp + 2Vg- Vp, 
we obtain from ( ^.6| ) the equation 

A(g/) = pAg + gAp + 2Vg • Vp + /^{qg). 
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Let's consider the degree, as a polynomial in Xr+i, . . . , a;„, of each term on the right 
side of this equation. Because Aq is a positive constant, this degree is k for the 
first term. Because q is independent of the variables Xr+i, . . . the second and 
third terms have degree (as a function of Xr+i, ■ ■ ■ ,Xn) less than k. Because the 
degree of g (as a function of x^+i, . . . , is less than k and q is independent of 
Xr+i, . ■ . ,Xm the degree of the fourth term is less than k. 

Thus the only part of the right side of the equation above with degree k (as 
a function of Xr+i, ■ ■ ■ ,Xn) is the first term, pAq. The left side of the equation 
is 0, so pAq — 0. Hence p — 0. But p was the part of / of highest degree in 
Xr+i, . . . ,Xn- Hence / is independent of Xr+i, . . . ,Xn- Thus we can think of the 



equation A{qf ) = as taking place in R*". Lemma 2.3 now implies that / = 0, as 



desired. □ 



We now apply Theorem 2.5 to obtain the general decomposition theorem. 



Theorem 2.8. Suppose p G Vm- Let q be a nonhyperbolic quadratic. Then there 
exists a unique harmonic polynomial h G Vm such that 

p = h + qf 

for some f G Vm~2- 



Proof. Note that Ap G Vm-2- Thus by Theorem 2.5, there exists / G Pm-2 such 
that A{qf) = Ap. Let h ~ p — qf . Then h is harmonic polynomial in Vm and 
p =: h + qf, as desired. 

To prove the uniqueness part of this theorem, suppose also that ft, is a harmonic 
polynomial in Vm and that p = h + qf for some / G Vm-2- Then 

h-h^q{f-f). 

The left side of the equation above is harmonic, and hence A(g(/ — /)) = 0. 



Theorem 2.5 now implies that / — / = 0, which implies that /i = ft, as desired. □ 



If p G Vm and g is a nonhyperbolic quadratic, we can consider the following 
Dirichlet problem: find a harmonic polynomial ft G Vm such that ft equals p on 
the set {x G R" : q{x) = 0}. Clearly the ft produced by the theorem above solves 
this Dirichlet problem. Of course, {x G R" : q{x) — 0} could be the empty set 
or a single point. Because q{x) — > cx) as \x\ — > oo, the existence of a point in R" 
where q is negative is a convenient condition to ensure that {x G R" : q{x) = 0} is a 



nondegenerate quadratic surface. For example, using the notation of Definition 2.4, 
{x G R" : q{x) = 0} will be a nondegenerate quadratic surface if 

^< S 4 

or if Cj ^ for some j with bj =0. 

We now turn to the question of whether the polynomial h produced by The- 
orem |2.8| is the unique polynomial solution to the Dirichlet problem discussed in 
the paragraph above. The following lemma will help us answer this uniqueness 
question. If we were working in C" instead of R", then Hilbert's NuUstellensatz 
could be used to provide information about when a polynomial ft vanishing on the 
zero set of another polynomial g is a polynomial multiple of q. A theorem called 
the Real NuUstellensatz (see, for example, Chapter 3, Theorem 3.3) provides 
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some information about polynomials on R" vanishing on the zero set of another 
polynomial. However, we do not see how the Real Nullstellensatz can be used to 
prove the lemma below, so we have provided a proof without using such machinery. 

Lemma 2.9. Suppose q is a nonhyperholic quadratic that is negative at some point 
ofW^. If h is a polynomial on R" such that h{x) = whenever q{x) = 0, then h 
is a polynomial multiple of q. 

Proof. We will prove this lemma by induction on the dimension n. 

To get started, suppose n = 1 and that q{x) = h'^x^ + ex d, where 6^0, 
is a nonhyperholic quadratic that is negative at some point of R. Suppose ft is a 
polynomial on R such that h{x) = whenever q{x) = 0. Because q is negative at 
some point of R and ^ cx) as a; — > oo, we see that q has precisely two distinct 
zeros. The polynomial h vanishes on both these zeros, and thus /i is a polynomial 
multiple of the quadratic polynomial q, as desired. 

Now suppose that the lemma holds in dimension n—1. Let g be a nonhyperholic 
quadratic that is negative at some point of R". Relabelling coordinates, if necessary, 
we can assume that 

n n 

where bj ^ for some j E {1, . . . , n — 1}. Let y denote a typical point of R"^^ and 
let z denote a typical point of R; thus (j/, z) denotes a typical point of R". 

Suppose /i is a polynomial on R" of degree m such that h{x) = whenever 
q{x) = 0. We need to show that /i is a polynomial multiple of q. 

For z gR such that q{y, z) is negative for some y e R"~^, define a nonhyperholic 
quadratic on R"~^ by 

<iz{y] = q{y,z), 

and define a polynomial hz on R"~^ of degree at most m by 

hz{y) = h{y,z). 

Then hz{y) = whenever qz{y) = 0, and thus by our induction hypothesis there is 
a polynomial fz on R"~^ such that 

(2.10) hz{y) ^ fz{y)qz{y) 

for all y € R"~^. Clearly fz has degree at most m — 2. 

For y e R"~^ such that q{y, z) is negative for some 2; G R, define a polynomial 
qy on R by 

qy{z) = q{y,z), 
and define a polynomial on R of degree at most m by 

hy{z) = h{y,z). 

Then hy{z) = whenever q^{z) = 0. We claim that there is a polynomial on R 
such that 

(2.11) hy{z) = gy{z)qy{z) 

for all 2: G R. If ^ 0, then q^ is a nonhyperholic quadratic on R and the claim 
follows from the dimension 1 case that was proved at the beginning of this proof. If 
bn = but Cn ^ 0, then q^ is a polynomial on R of degree 1 and the claim follows 
easily. Finally, if 6„ = and c„ = 0, then q^ is a negative constant, in which case 
the claim is trivially true. In any case, we see that g^ has degree at most m. 
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Let n = {{y,z) eK^ : z) < 0}. Combining ( |2.10D and ( |2.1l| ), we see that 

Hy, z) = fz{y)q{y, z) = gy{z)q{y, z) 
for all {y, z) e i7. Thus we can define a function p on by 

p{y,z) = f,Xy) = ^^(z). 

Because 

for all (y, z) e f2, we see that p is real-analytic on fi. 

Suppose a = (ai, . . . , a„) is a multi-index of order greater than 2m. Then either 

q;i + • ■ • + > m — 1 or a„ > m -|- 1. 

This implies that D°'p{y, 2;) = for all {y, z) G ft. Because all sufficiently high- 
order partial derivatives of p equal on the open set fl, we conclude that p is a 
polynomial on rt. Hence we can think of p as a polynomial defined everywhere 
on R". 

Finally, because h{y, z) — p{y, z)q{y, z) for all {y, z) G O, and because polynomi- 
als that agree on a nonempty open subset of R" must agree everywhere, we have 
h = pq. Thus h is a, polynomial multiple of g, as desired. □ 



Now we can combine the previous lemma and Theorem 2.5 to prove the desired 
uniqueness result. Of course in the ellipsoidal case (where each bj ^ 0, in the 



notation of Definition 2.4) uniqueness follows easily from the boundedness of the 
surface, but we want to consider also elliptic cylinders and paraboloids. Note that 
the uniqueness result in the theorem below fails on some nondegenerate quadratic 
surfaces, so the hypothesis that q is nonhyperbolic cannot be deleted. For example, 
on the quadratic surface defined by {x G R" : — 3x2 — 1 = 0}, any solution to 
a Dirichlet problem can be used to produce another solution by adding to it the 
harmonic polynomial [x^ — 3a;| — l)xi, which vanishes on the quadratic surface in 
question. 

To obtain uniqueness results on half-spaces, even in the class of polynomial 
solutions, a growth condition on the solutions is needed (see j7|). However, the 
theorem below shows that we have unique polynomial solutions on our quadratic 
surfaces without the requirement of a growth condition. 

Theorem 2.12. Suppose q is a nonhyperbolic quadratic that is negative at some 
point o/R". If p G Vrn, then there is a unique harmonic polynomial h that equals 
p on {x G R" ; q{x) = 0}. Furthermore, 

h=p-qf 

for some f G Vm-2- 



Proof. Take h and / as in Theorem 2.8. Then ft, is a harmonic polynomial that 



equals p on {x G R" : q{x) = 0}; furthermore / G 'Pm-2 and h—p — qf. 

To prove uniqueness, suppose h is also a harmonic polynomial that equals p on 



{x G R" : qix) = 0}. Then h-h equals on {x G R" : q{x) = 0}. By Lemma 2.9 



ft — /i is a polynomial multiple of q. However, Theorem 2.5 implies that no nonzero 
polynomial multiple of q is harmonic. Thus ft — ft = and hence h = h, completing 
the proof of uniqueness. □ 
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3. Algorithm 

In this section, we will turn the results of the previous section into a compu- 
tationally useful algorithm for solving the Dirichlet problem. Every polynomial 
is a sum of homogeneous polynomials, and thus to solve Dirichlet problems with 
polynomial boundary functions it suffices to do so for homogeneous polynomials. 

Let Hm denote the vector space of polynomials on R" that are homogeneuos 
of degree m (the polynomial is homogeneous of every degree, so G Hm)- Sup- 
pose g is a nonhyperbolic quadratic and p € Ti.m+2- We use the decomposition in 



Theorem 2.8 to write p ^ h + qf , where /i is a harmonic polynomial in Vm+2 and 
/ G Vm- Breaking each polynomial into its homogeneous components, we obtain 

m+2 m 

P=^^] + (92 + qi+ qo) fj: 

3=0 j=0 

where each hj is harmonic. We can break the equation above into homogeneuous 
equations by degree to obtain the system 

P = h,n+2 + q2fm, 
— Qlfm, — Kn+1 + 92/771-1, 
^Qofni — qifm-1 = Kn + <Z2/m-2, 

(3.1) : 

-90/2 - q1.f1 = h2+ 92/0, 
-90/1 - qifo = hi, 
-qofo = ^0- 

Here p and q2,qi, qo are known and we need to compute hm+2 , ■ ■ ■ ,ho and fm , . . . , /o . 
Our plan of attack is to use the first equation to find fm, which will then give us 
both hjn+2 and the left side of the next equation. We can repeat the procedure with 
the second equation to find /m-i, which will give us both hm+i and the left side of 
the next equation. We continue this process until we have found hj for each j. Then 
h = hj is the harmonic function that agrees with p on {a; e R" : q{x) — 0}. 

Thus, we turn our attention to equations of the form 

(3.2) P = h + q2f 

where p is a known homogeneous polynomial, q2 is the highest degree part of a 
known nonhyperbolic quadratic, h is an unknown harmonic polynomial homoge- 
neous of degree deg(p), and / is an unknown polynomial homogeneous of degree 
deg(p) — 2. Note that all the equations in the system above are of this form (except 
the last two equations, which are trivial to solve for hi and Hq). Hence an algorithm 
for finding the solution to ( |3.2| ) will give us an algorithm for solving our Dirichlet 
problem. 

To solve equation (^), eliminate h by taking the Laplacian of both sides, getting 

(3.3) Ap = A{q2f). 

Now our problem has been reduced to finding / e Tim satisfying the equation 
above, where p G 'Hm+2 an d q2{ x) = h\x\ -I- • • ■ -I- are known, with at least one 
hj / 0. Note that Theorem implies that (O) has a unique solution / e Hm- 
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The algorithm wc introduce for finding / involves repeated differentiation of ( 3.2 ). 
Because / S Hm, we will know / once we know the constants D"f for every multi- 
index a of order to. We need one more piece of notation before finding the constants 
D"f. For 1 < j < n, let ej denote the multi-index whose j"^-coordinate equals 1 
and whose other coordinates equal 0. 

For a a multi-index of order to, apply the differential operator 13" to both sides 



of (3.3), getting 



(3.4) D^{Ap)^D^{A{q2f)) 
= A(i?"(g2/)) 



where the last equality comes from Proposition A. 2 in our Appendix and the explicit 
form of • Because D"f is a constant, the Laplacian of the first term in parentheses 
above equals 2||6|pD"/. Because D'^~'^^ f has degree 1, the Laplacian of the second 
term in parentheses above can be easily computed using the product formula for 



the Laplacian (2.7), and the last equation becomes 



(3.5) Z3"(Ap) - (2| 



4E< 



+2efc-2e. 



/)]■ 



Note that multi-index a + 2ek — 2ej in the equation above has order to. Thus as 



a ranges over all multi-indices of order to, ( |3.5[ ) gives us a system of linear equations 
in the unknowns {D^f}. This system of equations can be solved using Gaussian 
elimination, giving us / and thus solving our Dir ichlet problem. 

Any solution to the syste m o f equations ( |3.5|) in the unknowns {D"f} gives a 
function / G Tim satisfying (3^) for all multi- indices a of order m, which implies 
that / satisfies Ap ~ A{q2f). However, we already know from Theorem that 
there is a unique / G Tim satisfying Ap = A{q2f). Thus the system of equations 
( |3.5| ) in the unknowns {D"/} has a unique solution. 

The number of operations needed to compute the solution to a system of s linear 
equations in s unknowns is on the order of (2/3)s'^. Thus we can expect that solving 
the system (3.5) will take on the order of 



(3.6) 



2(m- 



1)!3 



3TO!-^(n — 1)! 



operations (the formula for the number of multi-indices of order m can be found, 
for example, on page 78 of 0). However, a careful look at the system (3^) shows 
that we can do much better. 

For a fixed multi-index a, the multi-indices that appear in (^^) all have the 
form a + 2ek ~ 2ej. This leads us to define an equivalence relation on the set of 
multi-indices of order to by declaring that two multi-indices a and (3 are equivalent 
if ai = Pi (mod 2) for each i = 1, . . . , n. This equivalence relation breaks the set of 



multi- indices of order m into equivalence classes, and the system of equations (3.5) 
breaks into corresponding systems of equations. Hence instead of solving one large 
system of equations, we can solve several smaller systems of equations, which leads 
to considerable computational savings, as we will soon see. 
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If in is larger than n, then we have 2"^^ equivalence classes, corresponding to 
a choice of even or odd entry in the first n ~ 1 coordinates of a multi-index (the 
parity of the last coordinate is forced by the condition that the coordinates add up 
to m). As TO increases, the number of elements in each of these 2"~^ equivalence 
classes divided by the total number of multi-indices of order to approaches 1/2"^^. 
Thus instead of solving one large system of {"^~^^~ ) equations in the same number 
of variables, we can solve 2"~^ systems of equations, each system containing on the 
order of (™^'^^) /2"~^ equations and variables. This computation requires on the 
order of 

.o7^ 1 2(m + n-l)!3 

^ ■ ' 3(2"-i)3TO!3(n- 1)!3 

operations. 



The ratio of (3.6) to (3.7) is 2^"^^. Thus breaking our system of equations into 
smaller systems of equations using our equivalence relation reduces the number of 
operations by a factor of 2^"^^. For example, if n = 6, then this technique should 
reduce computation time by a factor of over 1000. This savings is needed even for 



moderate size to and n, because the system of equations (3.5) grows large rapidly. 

Usually we can do even better than reducing computations by a factor of 2-^"^^. 
Suppose, for example, that p{x) = xf^xl- Then Ap — 42x1^X2 + 380x}^a;2. 



Note that the left side of (3.5) equals for all multi-indices of order 25 except 
(20, 5, 0, . . . , 0) and (18, 7, 0, ... , 0). In other words, in only one of our equivalence 
classes (the equivalence class consisting of those multi-indices of order 25 whose 
first coordinate is even, second coordinate is odd, and coordinates 3 through n are 
even) is the left side of ( ^.5D anything other than 0. When the left side of every 
equation in an equivalence class is 0, there is no need to perform Gaussian elimina- 
tion to solve the system of equations in that equivalence class, because obviously 



all the unknowns equal (recall that the system (3.5) has a unique solution). Thus 
in the example at hand, instead of solving 2""-'^ smaller systems of equations, we 
need only solve one smaller system of equations. 

As can be seen from the reasoning in the previous paragraph, if p is a monomial 
then our computation time is reduced, through the use of our equivalence classes, 
by another factor of 2"~^. Thus if p is a monomial, then computation time through 
the use of equivalence classes is reduced by a factor of 2'^"^^. For example, if n = 6, 
then this technique should reduce computation time by a factor of over 32,000. 

The full strength of the reduction discussed above by a factor of 2^"^^, as opposed 
to a reduction by a factor of 2^"^^, holds only for ellipsoids and elliptic cylinders. 
For ellipsoids and ell ipti c cylinders we can assume (perhaps after a tran slation) that 



each Cj in Definition |2.4| equals 0, which gives = in the system { p.l\) . However, if 



qi ^ then multiplication by qi , when solving the second and successive equations 



in the system (3.1), can lead to nonzero left sides in equations in the system (3.5) 
other than the equations in the equivalence class corresponding to the exponents 
in the monomial p. 

4. Examples 

The algorithm outlined in the previous section has been implemented by the 
authors in Mathematica and in MATLAB. The Mathematica version produces exact 
solutions, in considerably less time than we believed possible when we started this 
project. Even with simple boundary polynomials in low dimensions, solutions to 
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Dirichlct problems on nonhyperbolic quadratic surfaces tend to involve fractions 
with largo numerators and denominators, as we will see in the examples presented 
below. Our Mathematica software can use floating point arithmetic to produce 
only decimal approximations to the exact solutions, with even faster times than 
when working with exact rational arithmetic. Our MATLAB implementation of 
the algorithm works only in floating point arithmetic, again quickly producing 
decimal approximations to the exact solutions. 

Our Mathematica implementation of the algorithm is available on the first au- 
thor's web site; our MATLAB implementation of the algorithm is available on the 
third author's web site. The appropriate web addresses are listed at the end of 
this paper; within those web sites look for information about this paper to find 
the software. Although our software is available without charge, its use requires 
Mathematica or MATLAB. The examples presented below were generated by our 
Mathematica implementation of the algorithm. 

Our first example will be in dimension 3 with a boundary function of degree 7. 
Here, as elsewhere throughout this section, all fractions are given in reduced form. 

Example 4.1. Suppose p{xi,X2,x^) 

q{xi,X2,X2,) = 2x1 + 3x1 +4x3 - 1. 

Then the following function is harmonic on and agrees with p on the ellipsoid 
{x€B? : q{x) = 0}; 

43 2 „ 2 .2 97950 5 2524856930 2 

x^xi, + (2xi + 3.T0 + Axi - 1) Xo xix2 

1 2 -TV i-r 2-r 3 ''V 20144813 ^ 100139865423 ^ 

3423451 4 148091 3 2306686 2 3 701980831 



X1X2 — 7r7r7;7;7T7T7r77T X2 — ^^010 ^1^2 ^ X2 



60434439 ' 33379955141 ^ 20144813 ' " 500699327115 

32326712 2 3712712 2 2 53836 3 2 _ 236464 4\ 
^ 7703066571 ^'^^ ^ 60434439 ^^^'^^ + 20144813 ^^^^ ~ 60434439 ^'^V ' 



The function above obviously equals xfxg on the ellipsoid in question. Thus 
to verify that the function above is indeed the solution to our Dirichlet problem, 
it is only necessary to verify that the Laplacian of the function above equals 0; 
Mathematica can easily perform this calculation. 

Note that, as expected from our discussion in the previous section, every ex- 
ponent of each Xj in the solution above has the same parity as in the boundary 
function. Specifically, in the solution above the exponents of Xi, x^, and Xi are 
even and the exponent of X2 is odd, thus following the pattern of the boundary 
function x\x\. 

Our Mathematica implementation of the algorithm can handle quadratic surfaces 
defined with symbols as well as concrete numbers. This capability is illustrated in 
our second example, which takes place in dimension 4 with a boundary function of 
degree 7. 

Example 4.2. Suppose p{x\,X2, Xi,Xi) 



q{xi,X2, X3, X4) = cx\ + 3x2 + 4x3 + 5x4 — 1. 
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Then the following function is harmonic on lA^ and agrees with p on the ellipsoid 
{x e R-* : q{x) = 0}; 

3 2 



788400 + 367920c + 48712c2 + 2520c3 + 45c' 



(4(50 + 3c) (36 + bc)xix\xji 



.o.o.nn no. n2^ 3 (82800 + 21868c + 1764c2 + 450^) 

- 12(2190 + 281c + %c^)x\x^Xi - — ^ XiX^^Xa 

3(10 + cj 

- (50400 + 21118C+ 1845c2 + Abc^)xixlx:iXA + 5(46 + 3c)(36 + bc)xiXzxl 

We now return to dimension 3 but increase the degree of the boundary function 
to 10. 

Example 4.3. The harmonic polynomial that equals x^* on the ellipsoid 
{xeH^ : 2x1 + 3x1 + ixj -1 = 0} 

has value 

500945213823452554440546462385400584789 



397263369506735959801289842040922215251461 

at the origin. 

Of course, the value at the origin of the solution to this Dirichlet problem is 
simply the constant term in the polynomial that gives the solution. In the example 
above, we have given only the constant term because displaying the entire solution 
would require more than a page. 

The example above illustrates the difference between the sphere and ellipsoids. 
The solution to the Dirichlet problem on the unit sphere in R'^ with boundary 
function has value 1/11 at the origin, and the largest integer appearing in the 
numerator or denominator of the coefficients of any term of the solution is 46189. 
In contrast, if the ellipsoid {x G R'^ : 2x1 + + "^^3 — 1 = 0} replaces the unit 
sphere, then the value of the solution at the origin is given by the fraction in the 
example above, and the coefficients of the other terms of the solution have similarly 
huge numerators and denominators. 

We conclude this section by giving some data about the speed of our algorithm. 
The times given below are CPU times used by the Mathematica implementation of 
our software on a Windows desktop computer with a 1.7 GHz Intel Pentium 4 chip. 



Examples 4.1 and 4.2 each took about 0.04 seconds. To produce Example 4.3 
we found the solution to the Dirichlet problem with boundary function xY^ on the 
ellipsoid {a; S R'^ : 2x\ + 3x2 + ~ 1 = 0}; this took less than 0.2 seconds. On 
the same ellipsoid, solving the Dirichlet problem with boundary function xf^ takes 
55 seconds. Changing the boundary function to x\^ increases the solution time to 
10 minutes, 5 seconds. Finally, changing the boundary function to x^ (still on the 
same ellipsoid) increases the solution time to 2 hours, 29 minutes, 32 seconds. 

The long solution times reported in the paragraph above arise because of the 
huge overhead associated with exact manipulation of the gigantic numerators and 
denominators that appear in the solutions to these high-degree problems. For exam- 
ple, most of the numerators and denominators are larger than 10^°" for the Dirichlet 
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problem above with boundary value xf", rising to over 10^^*^ for the Dirichlet prob- 
lem above with boundary value xl^. 

However, our algorithm is fast even for high-degree boundary functions if we use 
Mathematical s floating point arithmetic instead of exact rational arithmetic. For 
example, if we ask the Mathematica implementation of our algorithm to use floating 
point arithmetic to solve the Dirichlet problem mentioned above with boundary 
function xf", then the computation time is reduced from almost 2.5 hours to find 
the exact solution to just 8.9 seconds to find a very good decimal approximation of 
the solution. 



Appendix 

This purpose of this appendix is to prove the differentiation formula given by 
Proposition A. 2, which was used in the derivation of our algorithm in Section ^ 

We begin with following lemma, which gives a formula for differentiating the 
product of a linear function and a polynomial. 

Lemma A.l. Suppose f is a polynomial on R" and a is a multi-index. Then 



for every g Vi. 



Proof. Fix g Vi. We will prove our desired result by induction on jaj. To get 
started, note that the desired result is obviously true when |q;| = 0. 

Now suppose that |q;| > and that the desired result holds for all niulti- indices 
of smaller order. Choose k such that ak > 0. Then 

D^{gf) = D''-'^-Dk{gf) 

- D'^-'- {gDkf + {Dkg)f) 

= D''-'^^{gDkf) + {Dkg){D"-^''f), 

where the last equation holds because Dkg is constant. Applying our induction 
hypothesis to evaluate D°'~^''{gDkf) now gives 

n 

D"{gf) = gDy + ^ (Z?,,g)(D"-^^ /), 

completing the proof. □ 

The next proposition, which gives a formula for differentiating the product of a 
quadratic polynomial and another polynomial, was used in deriving the system of 
equations (3.5). 

Proposition A. 2. Suppose f is a polynomial on R" and a is a multi-index. Then 

n 1 ^ 

D'^iqf) = qOy + J2a,{D,q){D--^^f) + - 5]a,(a, - l){D]q){D--'^^ f) 
for every nonhyperbolic quadratic q. 
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Proof. Fix a nonhypcrbolic quadratic q. We will prove our desired result by in- 
duction on I a I . To get started, note that the desired result is obviously true when 
|a| = 0. 

Now suppose that |a| > and that the desired result holds for all multi- indices 
of smaller order. Choose k such that > 0. Then 

D''{qf)=D''-''^Dk{qf) 

(A.3) = D'^-'" (qDkf) + D^--" {{Dkq)f) . 



Using Lemma A.l to evaluate the last term (applicable because Dkq G "Pi) gives 
(A.4) D'^-^>'{{Dkq)f) = iDkq)iD'^-'''f) + iak-l){Dlq){D'^-^'^''f), 

where we have used the fact that DjO^q = whenever j ^ k. Using our induction 



hypothesis to evaluate the first term on the right side of (A.3) and using (A.4) to 



evaluate the second term on the right side of ( [A.3| ) now gives 

n 1 

D'^iqf) = qD'^f + ^a,(i?,g)(7?"-«V) + 3 E "^("^ " ^){D]q){D'^-^'^^ f), 
j=i i=i 
completing the proof. □ 
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